#ifndef _POTENTIAL_ENERGY_H_
#define _POTENTIAL_ENERGY_H_

#include <algorithm>
#include <cmath>
#include <vector>

#include <knitro.h>

#ifdef _CILK_
#include <cilk/cilk.h>
#endif

#include "containers.h"


#include "glut_viewer/glut_viewer.h"

namespace Config_parameters
{
	typedef enum { DIRECT = 1, CG = 2 } Opti_alg;
	extern Opti_alg opti_algo;
	extern double gravity_scale;
}

namespace Packing_common
{
	/*
	** C style smart pointer
	*/
	template <class T>
	class c_auto_ptr
	{
	public:
		explicit c_auto_ptr(T *p = 0) { ptr = p; }
		~c_auto_ptr() { if (ptr) free(ptr); }

		c_auto_ptr& operator=(c_auto_ptr& other) { ptr = other.ptr; other.ptr = 0;return *this; }
		c_auto_ptr& operator=(T *p) { ptr = p; p = 0; return *this; }
		T& operator*() const { return *ptr; }
		T* operator->() const { return ptr; }
		operator T*() const { return ptr; }

	private:
		T *ptr;
	};
	/*
	** Parameter collection for Knitro optimization
	** In order to reduce code length and,
	** avoid redudant setup of Knitro context whenever an optimzation is executed
	*/
	struct Knitro_context
	{
		Knitro_context(int nb_vars, int nb_constraits, int nb_nonzero_jac, int nb_nonzero_hess) ;
		~Knitro_context();
		int init_problem();
		KTR_context *kc;
		int n; // number of variables
		int m; // number of constraints
		int nnzJ; // number of nonzero elements in Jacobian
		int nnzH; // number of nonzero elements in Hessian
		int objGoal;
		int objType;
		c_auto_ptr<int> cType; // constraint type
		c_auto_ptr<int> jacIndexVars; // 
		c_auto_ptr<int> jacIndexCons;
		c_auto_ptr<int> hessIndexRows;
		c_auto_ptr<int> hessIndexCols;
		double obj;
		c_auto_ptr<double> lambda;
		//double *xLoBnds;
		//double *xUpBnds;
		//double *cLoBnds;
		//double *cUpBnds;
		//double *xInitials;
		c_auto_ptr<double> xLoBnds;
		c_auto_ptr<double> xUpBnds;
		c_auto_ptr<double> cLoBnds;
		c_auto_ptr<double> cUpBnds;
		c_auto_ptr<double> xInitials;
	};

	class System_potential_energy
	{

	public:

		System_potential_energy(const Container_box& bndbox) : bounding_box(bndbox) { }
		virtual int optimize(double *trans, bool apply_constraints = false) = 0;

		// objective function
		virtual double system_energy_at(const double *const x) 
		{
			return gravity_energy_at(x) + molecular_energy_at(x);
		}
		virtual double gravity_energy_at(const double *const x) = 0;
		virtual double molecular_energy_at(const double *const x) = 0;

		virtual void energy_gradient_at(const double *const x, double *gradient) = 0;	

		virtual ~System_potential_energy() { }

	private:

		// constraint function
		virtual void constraint_at(const double *const x, double *constraint_values) = 0;
		virtual void constraint_grad_at(const double *const x, double *constraint_grad) = 0;

		// gradient of the objective function
		//virtual void energy_gradient_at(const double *const x, double *gradient) = 0;	

	protected:
		static int callback(const int evalRequestCode, const int n, const int m, const int nnzJ, 
							const int nnzH,	const double * const x,	const double * const lambda, 
							double * const obj, double * const c, double * const objGrad, double * const jac,
							double * const hessian,	double * const hessVector, void * userParams);

	protected:
		const Container_box& bounding_box;
		bool nonoverlap_constraint;

		static double epsilon;
		static double g; // for gravity potential
		static double f; // for adjusting minimal energy (r_m)
		//static double t; // adjustor of barrier function
	};


}
#endif